The General Linear Model

Extending Effects

Table of Contents


You have already seen the basics of the general linear model (linear regression, t-tests, and ANOVAs). Now, we are going to get into some expanded topics related to the general linear model: centering, interactions, and mediation.

Main Ideas

Centering – Makes interpretable intercepts

Interactions – Adds context to models

Mediation – Adds paths to models

Centering

Think back to what the intercept means in the context of a linear regression model: it is the value of the DV/outcome when everything else is at 0. There are times where that makes sense (you can have 0 dollars, but I don’t recommend it).

What if we want to have a meaningful intercept?


library(data.table)

crimeScore <- fread("http://nd.edu/~sberry5/data/crimeScore.csv")

uncentered <- lm(SSL_SCORE ~ NARCOTICS_ARR_CNT, data = crimeScore)

summary(uncentered)

Call:
lm(formula = SSL_SCORE ~ NARCOTICS_ARR_CNT, data = crimeScore)

Residuals:
     Min       1Q   Median       3Q      Max 
-255.974  -44.233    8.767   42.285  223.767 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       275.4925     0.3175  867.70  < 2e-16 ***
NARCOTICS_ARR_CNT   0.7408     0.1204    6.15 7.76e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 60.39 on 92805 degrees of freedom
  (305877 observations deleted due to missingness)
Multiple R-squared:  0.0004074, Adjusted R-squared:  0.0003967 
F-statistic: 37.83 on 1 and 92805 DF,  p-value: 7.758e-10

plot(uncentered$fitted.values, uncentered$model$SSL_SCORE)


plot(uncentered$model$SSL_SCORE, uncentered$model$NARCOTICS_ARR_CNT)

Now we can do some very simple centering: we can just subtract the mean of the item from every observation of that item.


crimeScore[, narcCenter := NARCOTICS_ARR_CNT - mean(NARCOTICS_ARR_CNT, 
                                                    na.rm = TRUE)]

# Notice how I didn't assign crimeScore back to crimeScore?
# The skull operator handled it for us!

centeredMod <- lm(SSL_SCORE ~ narcCenter, 
                  data = crimeScore)

summary(centeredMod)

Call:
lm(formula = SSL_SCORE ~ narcCenter, data = crimeScore)

Residuals:
     Min       1Q   Median       3Q      Max 
-255.974  -44.233    8.767   42.285  223.767 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 277.0178     0.1982 1397.34  < 2e-16 ***
narcCenter    0.7408     0.1204    6.15 7.76e-10 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 60.39 on 92805 degrees of freedom
  (305877 observations deleted due to missingness)
Multiple R-squared:  0.0004074, Adjusted R-squared:  0.0003967 
F-statistic: 37.83 on 1 and 92805 DF,  p-value: 7.758e-10

plot(centeredMod$fitted.values, centeredMod$model$SSL_SCORE)


plot(centeredMod$model$SSL_SCORE, centeredMod$model$narcCenter)

We can see that nothing really changes except for the intercept’s coefficient. Let’s think about what is being said here: with narcotics arrest count being a 0, we would expect the SSL score to be 277 on average. Since we centered that narcotics variable, though, the 0 is actually the mean of that variable (2.06).

Standardizing

We will get back to more interesting data in a second, but let’s look at some mtcars data first (we can’t possibly learn anything new about it):


cor(mtcars[, c("hp", "disp", "wt")])

            hp      disp        wt
hp   1.0000000 0.7909486 0.6587479
disp 0.7909486 1.0000000 0.8879799
wt   0.6587479 0.8879799 1.0000000

That is a pretty strong correlation. Let’s see how they behave in a model:


testMod <- lm(mpg ~ hp + disp + wt, data = mtcars)

Since those variable had pretty strong correlations, let’s look at the variance inflation factor for each:

\[\frac{1}{1-R^2_i}\]


summary(lm(hp ~ disp + wt, data = mtcars))

Call:
lm(formula = hp ~ disp + wt, data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-53.830 -33.345  -3.866  15.445 155.541 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  68.8443    31.8023   2.165 0.038776 *  
disp          0.5388     0.1350   3.990 0.000411 ***
wt          -14.4453    17.1038  -0.845 0.405268    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 42.85 on 29 degrees of freedom
Multiple R-squared:  0.6346,    Adjusted R-squared:  0.6094 
F-statistic: 25.18 on 2 and 29 DF,  p-value: 4.575e-07

We can plug our \(R^2\) value in:


1 / (1 - .6346)

[1] 2.736727

VIF starts at 1 and goes up from there. Anything over 5 would indicate worrisome multicolinearity.


car::vif(testMod)

      hp     disp       wt 
2.736633 7.324517 4.844618 

When modeling these main effects, nothing is going to change this problem. If we are going to be using our variable as interactions, we might want to standardize these variables – we would center our items and then divide by the standard deviation of the item:


mtcarsDT <- as.data.table(mtcars)

scaleVars <- c("hp", "disp", "wt")

mtcarsDT[, 
         (scaleVars) := lapply(.SD, function(x) scale(x)), 
         .SDcols = scaleVars]

summary(mtcarsDT)

      mpg             cyl             disp               hp         
 Min.   :10.40   Min.   :4.000   Min.   :-1.2879   Min.   :-1.3810  
 1st Qu.:15.43   1st Qu.:4.000   1st Qu.:-0.8867   1st Qu.:-0.7320  
 Median :19.20   Median :6.000   Median :-0.2777   Median :-0.3455  
 Mean   :20.09   Mean   :6.188   Mean   : 0.0000   Mean   : 0.0000  
 3rd Qu.:22.80   3rd Qu.:8.000   3rd Qu.: 0.7688   3rd Qu.: 0.4859  
 Max.   :33.90   Max.   :8.000   Max.   : 1.9468   Max.   : 2.7466  
      drat             wt               qsec             vs        
 Min.   :2.760   Min.   :-1.7418   Min.   :14.50   Min.   :0.0000  
 1st Qu.:3.080   1st Qu.:-0.6500   1st Qu.:16.89   1st Qu.:0.0000  
 Median :3.695   Median : 0.1101   Median :17.71   Median :0.0000  
 Mean   :3.597   Mean   : 0.0000   Mean   :17.85   Mean   :0.4375  
 3rd Qu.:3.920   3rd Qu.: 0.4014   3rd Qu.:18.90   3rd Qu.:1.0000  
 Max.   :4.930   Max.   : 2.2553   Max.   :22.90   Max.   :1.0000  
       am              gear            carb      
 Min.   :0.0000   Min.   :3.000   Min.   :1.000  
 1st Qu.:0.0000   1st Qu.:3.000   1st Qu.:2.000  
 Median :0.0000   Median :4.000   Median :2.000  
 Mean   :0.4062   Mean   :3.688   Mean   :2.812  
 3rd Qu.:1.0000   3rd Qu.:4.000   3rd Qu.:4.000  
 Max.   :1.0000   Max.   :5.000   Max.   :8.000  

scaledMod <- lm(mpg ~ hp + disp + wt, data = mtcarsDT)

summary(scaledMod)

Call:
lm(formula = mpg ~ hp + disp + wt, data = mtcarsDT)

Residuals:
   Min     1Q Median     3Q    Max 
-3.891 -1.640 -0.172  1.061  5.861 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  20.0906     0.4665  43.067  < 2e-16 ***
hp           -2.1362     0.7841  -2.724  0.01097 *  
disp         -0.1161     1.2827  -0.091  0.92851    
wt           -3.7190     1.0432  -3.565  0.00133 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.639 on 28 degrees of freedom
Multiple R-squared:  0.8268,    Adjusted R-squared:  0.8083 
F-statistic: 44.57 on 3 and 28 DF,  p-value: 8.65e-11

Interactions

Sticking with some boring data for just a little while longer:


library(ggplot2)

ggplot(mtcars, aes(mpg, hp)) + 
  geom_point() +
  geom_smooth(method = "lm") +
  theme_minimal()

We can clearly see the effect of horsepower on mpg.

Let’s add some additional context to this visualization:


ggplot(mtcars, aes(mpg, hp, color = as.factor(cyl))) + 
  geom_point() +
  geom_smooth(method = "lm", se = FALSE) +
  theme_minimal()

We have already been doing stuff like this, but not really thinking about the model that could test it.

Let’s start by looking at a multiple regression with 2 predictor variables:


twoVars <- lm(SSL_SCORE ~ WEAPONS_ARR_CNT + NARCOTICS_ARR_CNT, 
             data = crimeScore)

summary(twoVars)

Call:
lm(formula = SSL_SCORE ~ WEAPONS_ARR_CNT + NARCOTICS_ARR_CNT, 
    data = crimeScore)

Residuals:
     Min       1Q   Median       3Q      Max 
-194.221  -28.029   -2.221   26.550  187.586 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       291.5470     1.5197 191.850   <2e-16 ***
WEAPONS_ARR_CNT    18.0598     1.0118  17.849   <2e-16 ***
NARCOTICS_ARR_CNT   2.8072     0.2933   9.573   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 49.69 on 9037 degrees of freedom
  (389644 observations deleted due to missingness)
Multiple R-squared:  0.04265,   Adjusted R-squared:  0.04243 
F-statistic: 201.3 on 2 and 9037 DF,  p-value: < 2.2e-16

Let’s explore interactions (moderation to some). A key idea here is that interactions change the nature of our model. Instead of supposing that the predictor variables act in isolation on the DV, the interaction is essentially providing expanded context for our model. We are essentially saying that the values of one variable will have an effect with values of another variable on the DV. Here is what this looks like in words:

“What is the effect of X on Y?”

“It depends on Z’s value.”


intMod <- lm(SSL_SCORE ~ WEAPONS_ARR_CNT * NARCOTICS_ARR_CNT, data = crimeScore)

summary(intMod)

Call:
lm(formula = SSL_SCORE ~ WEAPONS_ARR_CNT * NARCOTICS_ARR_CNT, 
    data = crimeScore)

Residuals:
     Min       1Q   Median       3Q      Max 
-194.231  -27.995   -2.231   26.532  187.532 

Coefficients:
                                  Estimate Std. Error t value
(Intercept)                       292.1106     2.2340 130.754
WEAPONS_ARR_CNT                    17.5941     1.6896  10.413
NARCOTICS_ARR_CNT                   2.5461     0.8134   3.130
WEAPONS_ARR_CNT:NARCOTICS_ARR_CNT   0.2173     0.6312   0.344
                                  Pr(>|t|)    
(Intercept)                        < 2e-16 ***
WEAPONS_ARR_CNT                    < 2e-16 ***
NARCOTICS_ARR_CNT                  0.00175 ** 
WEAPONS_ARR_CNT:NARCOTICS_ARR_CNT  0.73069    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 49.69 on 9036 degrees of freedom
  (389644 observations deleted due to missingness)
Multiple R-squared:  0.04266,   Adjusted R-squared:  0.04234 
F-statistic: 134.2 on 3 and 9036 DF,  p-value: < 2.2e-16

This is the model that we have estimated:

\[score = b_0 + b_1(weapons) + b_2(narcotics) + b_3(weapons*narcotics)\]

The interpretation of our main effects don’t really change.

Our interaction terms (\(b_3\)) is providing the amount of change in the slope of the regression of score on weapons when narcotics changes by one unit. So as narcotics arrests increase, we see an increase in the effect of weapons arrests on score (but only a tiny one).

To predict what the score value would be for certain values of weapons arrests, we could reformulate our model as:

\[score = b_0 + (b_1 + b3*narcotics)weapons + b_2(narcotics)\]

Before we look at those interactions, let’s compare how our models fit:


anova(twoVars, intMod)

Analysis of Variance Table

Model 1: SSL_SCORE ~ WEAPONS_ARR_CNT + NARCOTICS_ARR_CNT
Model 2: SSL_SCORE ~ WEAPONS_ARR_CNT * NARCOTICS_ARR_CNT
  Res.Df      RSS Df Sum of Sq      F Pr(>F)
1   9037 22311153                           
2   9036 22310861  1    292.56 0.1185 0.7307

We can compare some information about the residuals to compute an F statistic and the accompanying p-value – these models are no different from each other.


stargazer::stargazer(twoVars, intMod, type = "html", header = FALSE)
Dependent variable:
SSL_SCORE
(1) (2)
WEAPONS_ARR_CNT 18.060*** 17.594***
(1.012) (1.690)
NARCOTICS_ARR_CNT 2.807*** 2.546***
(0.293) (0.813)
WEAPONS_ARR_CNT:NARCOTICS_ARR_CNT 0.217
(0.631)
Constant 291.547*** 292.111***
(1.520) (2.234)
Observations 9,040 9,040
R2 0.043 0.043
Adjusted R2 0.042 0.042
Residual Std. Error 49.688 (df = 9037) 49.690 (df = 9036)
F Statistic 201.280*** (df = 2; 9037) 134.213*** (df = 3; 9036)
Note: p<0.1; p<0.05; p<0.01

library(effects)

modEffects <- effect("WEAPONS_ARR_CNT*NARCOTICS_ARR_CNT", intMod)

plot(modEffects)

This is offering the relationship between weapons arrests and score at various levels of narcotics arrests. What we are really looking for are slopes that differ from each other. Knowing that we didn’t find a significant interaction in our model, we probably shouldn’t be surprised that we don’t see any slope difference here.

The five levels that get plotted over is what would lead us towards something like a “floodlight analysis”. Below is more of a “spotlight analysis”:


library(interactions)

interact_plot(intMod, pred = WEAPONS_ARR_CNT, modx = NARCOTICS_ARR_CNT)


crimeScoreGender <- 
  crimeScore[SEX_CODE_CD != "X", .(SSL_SCORE, SEX_CODE_CD, WEAPONS_ARR_CNT)
  ][, SEX_CODE_CD := as.factor(SEX_CODE_CD)]

# Below is the dplyr equivalent to what we just used.
# Try it and see how slow it is...smh.

# crimeScoreGender = crimeScore %>% 
#   filter(SEX_CODE_CD != "X") %>%
#   dplyr::select(SSL_SCORE, SEX_CODE_CD, WEAPONS_ARR_CNT) %>% 
#   mutate(SEX_CODE_CD = as.factor(SEX_CODE_CD))

intMod2 <- lm(SSL_SCORE ~ WEAPONS_ARR_CNT * SEX_CODE_CD, data = crimeScoreGender)

summary(intMod2)

Call:
lm(formula = SSL_SCORE ~ WEAPONS_ARR_CNT * SEX_CODE_CD, data = crimeScoreGender)

Residuals:
     Min       1Q   Median       3Q      Max 
-264.092  -30.093   -0.093   29.908  182.908 

Coefficients:
                             Estimate Std. Error t value Pr(>|t|)    
(Intercept)                   316.467      7.939  39.862  < 2e-16 ***
WEAPONS_ARR_CNT                -2.251      7.202  -0.313  0.75465    
SEX_CODE_CDM                  -16.376      8.005  -2.046  0.04079 *  
WEAPONS_ARR_CNT:SEX_CODE_CDM   19.252      7.245   2.657  0.00788 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 53.03 on 19141 degrees of freedom
  (379482 observations deleted due to missingness)
Multiple R-squared:  0.02465,   Adjusted R-squared:  0.02449 
F-statistic: 161.2 on 3 and 19141 DF,  p-value: < 2.2e-16

Compared to women, men’s score increases by 19.25 on average for each weapons arrest.

Sometimes it helps to see what is going on with a plot:


modEffects <- effect("WEAPONS_ARR_CNT*SEX_CODE_CD", intMod2)

plot(modEffects)

Here is another way (and probably easier) to visualize this effect:


interact_plot(intMod2, pred = WEAPONS_ARR_CNT, modx = SEX_CODE_CD)

The visualization is telling, but the simple slope test is what will tell us which part of the model is significant:


sim_slopes(intMod2, pred = WEAPONS_ARR_CNT, modx = SEX_CODE_CD)

SIMPLE SLOPES ANALYSIS 

Slope of WEAPONS_ARR_CNT when SEX_CODE_CD = M: 

   Est.   S.E.   t val.      p
------- ------ -------- ------
  17.00   0.78    21.76   0.00

Slope of WEAPONS_ARR_CNT when SEX_CODE_CD = F: 

   Est.   S.E.   t val.      p
------- ------ -------- ------
  -2.25   7.20    -0.31   0.75

We see that the simple slope for F is not significant, whereas the simple slope for M is significant. This offers the “driver” of the relationship.

And another way:


sjPlot::plot_model(intMod2, type = "int", 
                    title = "") +
  ggplot2::theme_minimal()

Those plots offered the same information, just slightly different looks – pick whichever one you like most!

For the sake of it, let’s look at one more continuous by continuous interaction:


domNarc <- lm(SSL_SCORE ~ NARCOTICS_ARR_CNT * DOMESTIC_ARR_CNT, 
              data = crimeScore)

summary(domNarc)

Call:
lm(formula = SSL_SCORE ~ NARCOTICS_ARR_CNT * DOMESTIC_ARR_CNT, 
    data = crimeScore)

Residuals:
     Min       1Q   Median       3Q      Max 
-194.551  -30.857    6.657   32.449  224.421 

Coefficients:
                                   Estimate Std. Error t value
(Intercept)                        276.7127     1.0665 259.462
NARCOTICS_ARR_CNT                    2.8284     0.4019   7.038
DOMESTIC_ARR_CNT                    -1.4681     0.5447  -2.695
NARCOTICS_ARR_CNT:DOMESTIC_ARR_CNT   0.4775     0.2074   2.302
                                   Pr(>|t|)    
(Intercept)                         < 2e-16 ***
NARCOTICS_ARR_CNT                  2.02e-12 ***
DOMESTIC_ARR_CNT                    0.00704 ** 
NARCOTICS_ARR_CNT:DOMESTIC_ARR_CNT  0.02132 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 52.34 on 19090 degrees of freedom
  (379590 observations deleted due to missingness)
Multiple R-squared:  0.01304,   Adjusted R-squared:  0.01288 
F-statistic: 84.06 on 3 and 19090 DF,  p-value: < 2.2e-16

Let’s look at these interactions:


probe_interaction(domNarc, pred = NARCOTICS_ARR_CNT, modx = DOMESTIC_ARR_CNT)

JOHNSON-NEYMAN INTERVAL 

When DOMESTIC_ARR_CNT is OUTSIDE the interval [-49.01, -2.40],
the slope of NARCOTICS_ARR_CNT is p < .05.

Note: The range of observed values of DOMESTIC_ARR_CNT is [1.00,
24.00]

SIMPLE SLOPES ANALYSIS 

Slope of NARCOTICS_ARR_CNT when DOMESTIC_ARR_CNT = 0.47 (- 1 SD): 

  Est.   S.E.   t val.      p
------ ------ -------- ------
  3.05   0.33     9.33   0.00

Slope of NARCOTICS_ARR_CNT when DOMESTIC_ARR_CNT = 1.60 (Mean): 

  Est.   S.E.   t val.      p
------ ------ -------- ------
  3.59   0.23    15.67   0.00

Slope of NARCOTICS_ARR_CNT when DOMESTIC_ARR_CNT = 2.74 (+ 1 SD): 

  Est.   S.E.   t val.      p
------ ------ -------- ------
  4.14   0.33    12.52   0.00


sjPlot::plot_model(domNarc, type = "int", 
                    title = "") + 
  theme_minimal()

And also some categorical interactions:


crimeScore2 <- crimeScore[AGE_CURR != "" & SEX_CODE_CD != "X", 
           AGE_CURR_Rec := .(relevel(as.factor(AGE_CURR), ref = "less than 20"))]

# crimeScore2 = crimeScore %>% 
#   filter(AGE_CURR != "" & SEX_CODE_CD != "X") %>% 
#   mutate(AGE_CURR = relevel(as.factor(.$AGE_CURR), ref = "less than 20"))

sexMod <- lm(SSL_SCORE ~ SEX_CODE_CD, data = crimeScore2)

summary(sexMod)

Call:
lm(formula = SSL_SCORE ~ SEX_CODE_CD, data = crimeScore2)

Residuals:
    Min      1Q  Median      3Q     Max 
-273.46  -37.69   10.31   42.31  221.31 

Coefficients:
             Estimate Std. Error  t value Pr(>|t|)    
(Intercept)  283.4600     0.1868 1517.716   <2e-16 ***
SEX_CODE_CDM  -4.7709     0.2145  -22.246   <2e-16 ***
SEX_CODE_CDX -17.2845     7.6793   -2.251   0.0244 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 57.96 on 398681 degrees of freedom
Multiple R-squared:  0.001248,  Adjusted R-squared:  0.001243 
F-statistic:   249 on 2 and 398681 DF,  p-value: < 2.2e-16

ageMod <- lm(SSL_SCORE ~ AGE_CURR_Rec, data = crimeScore2)

summary(ageMod)

Call:
lm(formula = SSL_SCORE ~ AGE_CURR_Rec, data = crimeScore2)

Residuals:
     Min       1Q   Median       3Q      Max 
-256.449  -12.741   -1.449    9.567  200.551 

Coefficients:
                    Estimate Std. Error t value Pr(>|t|)    
(Intercept)        356.59475    0.09531  3741.5   <2e-16 ***
AGE_CURR_Rec20-30  -35.73399    0.10735  -332.9   <2e-16 ***
AGE_CURR_Rec30-40  -79.14579    0.11226  -705.0   <2e-16 ***
AGE_CURR_Rec40-50 -123.27298    0.12080 -1020.5   <2e-16 ***
AGE_CURR_Rec50-60 -166.16162    0.13211 -1257.7   <2e-16 ***
AGE_CURR_Rec60-70 -207.85339    0.19513 -1065.2   <2e-16 ***
AGE_CURR_Rec70-80 -254.62249    0.47812  -532.5   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 18.66 on 398379 degrees of freedom
  (298 observations deleted due to missingness)
Multiple R-squared:  0.8958,    Adjusted R-squared:  0.8958 
F-statistic: 5.708e+05 on 6 and 398379 DF,  p-value: < 2.2e-16

sexAge <- lm(SSL_SCORE ~ SEX_CODE_CD * AGE_CURR_Rec, data = crimeScore2)

summary(sexAge)

Call:
lm(formula = SSL_SCORE ~ SEX_CODE_CD * AGE_CURR_Rec, data = crimeScore2)

Residuals:
     Min       1Q   Median       3Q      Max 
-256.777  -12.739   -1.384    9.667  200.223 

Coefficients:
                                Estimate Std. Error  t value Pr(>|t|)
(Intercept)                     349.1487     0.1765 1978.133  < 2e-16
SEX_CODE_CDM                     10.4675     0.2093   50.018  < 2e-16
AGE_CURR_Rec20-30               -31.4097     0.2013 -156.055  < 2e-16
AGE_CURR_Rec30-40               -72.7648     0.2144 -339.444  < 2e-16
AGE_CURR_Rec40-50              -116.6520     0.2343 -497.896  < 2e-16
AGE_CURR_Rec50-60              -158.3222     0.2685 -589.698  < 2e-16
AGE_CURR_Rec60-70              -200.0051     0.4588 -435.955  < 2e-16
AGE_CURR_Rec70-80              -246.7251     1.2394 -199.063  < 2e-16
SEX_CODE_CDM:AGE_CURR_Rec20-30   -6.2596     0.2375  -26.356  < 2e-16
SEX_CODE_CDM:AGE_CURR_Rec30-40   -9.0747     0.2513  -36.112  < 2e-16
SEX_CODE_CDM:AGE_CURR_Rec40-50   -9.3964     0.2732  -34.400  < 2e-16
SEX_CODE_CDM:AGE_CURR_Rec50-60  -10.9607     0.3084  -35.540  < 2e-16
SEX_CODE_CDM:AGE_CURR_Rec60-70  -10.9464     0.5072  -21.581  < 2e-16
SEX_CODE_CDM:AGE_CURR_Rec70-80  -10.9950     1.3427   -8.189 2.65e-16
                                  
(Intercept)                    ***
SEX_CODE_CDM                   ***
AGE_CURR_Rec20-30              ***
AGE_CURR_Rec30-40              ***
AGE_CURR_Rec40-50              ***
AGE_CURR_Rec50-60              ***
AGE_CURR_Rec60-70              ***
AGE_CURR_Rec70-80              ***
SEX_CODE_CDM:AGE_CURR_Rec20-30 ***
SEX_CODE_CDM:AGE_CURR_Rec30-40 ***
SEX_CODE_CDM:AGE_CURR_Rec40-50 ***
SEX_CODE_CDM:AGE_CURR_Rec50-60 ***
SEX_CODE_CDM:AGE_CURR_Rec60-70 ***
SEX_CODE_CDM:AGE_CURR_Rec70-80 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 18.56 on 398372 degrees of freedom
  (298 observations deleted due to missingness)
Multiple R-squared:  0.8968,    Adjusted R-squared:  0.8968 
F-statistic: 2.664e+05 on 13 and 398372 DF,  p-value: < 2.2e-16

cat_plot(sexAge, pred = SEX_CODE_CD, modx = AGE_CURR_Rec, 
         geom = "line")

Mediation

As noted earlier, interactions are often referred to as moderators. Just to reiterate, they are testing if the relationship between the X & Y change when introducing a Z. A slightly different model details mediation – the relationship between X and Y exists because M directs the relationship.

We can refer to these paths as follows:

X -> M = a

M -> Y = b

X -> Y = c`

Let’s think how this might look conceptually:

Age -> Salary -> Job Satisfaction

We are saying that the effect of age on job satisfaction goes through salary. In other words, age doesn’t increase job satisfaction – age increases salary, which in turn increases job satisfaction.

To test this mediation model, we need to construct 4 separate regression models:

Model 1: Test path c (X predicts Y)

Model 2: Test path a (X predicts M)

Model 3: Test path b (M predicts Y)

If all of these yield significant results, you can go to:

Model 4: Use X and M as predictor to Y

If the effect of M on Y is significant in Model 4, then there is mediation: full mediation if X is not significant or partial mediation if X is significant. As awesome as a fully-mediated relationship sounds, it is not something that you are going to encounter all that often; instead, partial mediation is where things tend to fall.

Let’s try it out with our crime data:

Just to put some words behind this – we are proposing that narcotics arrests increase weapons arrests, which in turn leads to increased crime scores.


mediationData <- crimeScore[, 
                            .(SSL_SCORE, NARCOTICS_ARR_CNT, WEAPONS_ARR_CNT)]

mediationData <- na.omit(mediationData)

pathC <- lm(SSL_SCORE ~ NARCOTICS_ARR_CNT, 
            data = mediationData)

summary(pathC)

Call:
lm(formula = SSL_SCORE ~ NARCOTICS_ARR_CNT, data = mediationData)

Residuals:
    Min      1Q  Median      3Q     Max 
-198.43  -28.12   -2.43   26.57  183.26 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       314.0568     0.8627 364.060   <2e-16 ***
NARCOTICS_ARR_CNT   2.6866     0.2983   9.007   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 50.55 on 9038 degrees of freedom
Multiple R-squared:  0.008896,  Adjusted R-squared:  0.008787 
F-statistic: 81.13 on 1 and 9038 DF,  p-value: < 2.2e-16

Our pathC object (Model 1) is significant, so let’s go to the next step:


pathA <- lm(WEAPONS_ARR_CNT ~ NARCOTICS_ARR_CNT, 
            data = mediationData)

summary(pathA)

Call:
lm(formula = WEAPONS_ARR_CNT ~ NARCOTICS_ARR_CNT, data = mediationData)

Residuals:
    Min      1Q  Median      3Q     Max 
-0.2397 -0.2397 -0.2330 -0.2130  4.7670 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)        1.246405   0.008814 141.404   <2e-16 ***
NARCOTICS_ARR_CNT -0.006679   0.003048  -2.191   0.0285 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5165 on 9038 degrees of freedom
Multiple R-squared:  0.000531,  Adjusted R-squared:  0.0004204 
F-statistic: 4.802 on 1 and 9038 DF,  p-value: 0.02846

Model 2 – the effect of narcotics on weapons – is also significant. Things are looking promising so far.


pathB <- lm(SSL_SCORE ~ WEAPONS_ARR_CNT, 
            data = mediationData)

summary(pathB)

Call:
lm(formula = SSL_SCORE ~ WEAPONS_ARR_CNT, data = mediationData)

Residuals:
     Min       1Q   Median       3Q      Max 
-195.052  -28.052   -2.052   26.948  183.948 

Coefficients:
                Estimate Std. Error t value Pr(>|t|)    
(Intercept)      298.215      1.357  219.70   <2e-16 ***
WEAPONS_ARR_CNT   17.837      1.017   17.55   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 49.94 on 9038 degrees of freedom
Multiple R-squared:  0.03294,   Adjusted R-squared:  0.03283 
F-statistic: 307.8 on 1 and 9038 DF,  p-value: < 2.2e-16

Model 3 is also looking good! Weapons certainly has an effect on crime score.


pathFull <- lm(SSL_SCORE ~ NARCOTICS_ARR_CNT + WEAPONS_ARR_CNT, 
               data = mediationData)

summary(pathFull)

Call:
lm(formula = SSL_SCORE ~ NARCOTICS_ARR_CNT + WEAPONS_ARR_CNT, 
    data = mediationData)

Residuals:
     Min       1Q   Median       3Q      Max 
-194.221  -28.029   -2.221   26.550  187.586 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       291.5470     1.5197 191.850   <2e-16 ***
NARCOTICS_ARR_CNT   2.8072     0.2933   9.573   <2e-16 ***
WEAPONS_ARR_CNT    18.0598     1.0118  17.849   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 49.69 on 9037 degrees of freedom
Multiple R-squared:  0.04265,   Adjusted R-squared:  0.04243 
F-statistic: 201.3 on 2 and 9037 DF,  p-value: < 2.2e-16

Since all predictors continue to be significant, we would have partial mediation. It would be a weak moderation, though, because we did not reduced the magnitude of narcotics when the mediator was included in the model.

We can also test the indirect effect of X through M to Y by subtracting coefficients from models 1 and 4.


# This is Judd and Kenny's method:

pathC$coefficients["NARCOTICS_ARR_CNT"] - 
  pathFull$coefficients["NARCOTICS_ARR_CNT"]

NARCOTICS_ARR_CNT 
       -0.1206138 

Or by multiplying the full path moderator’s coefficient (model 4) from path a’s coefficient (model 2)


# This is Sobel's method:

pathFull$coefficients["WEAPONS_ARR_CNT"] * 
  pathA$coefficients["NARCOTICS_ARR_CNT"]

WEAPONS_ARR_CNT 
     -0.1206138 

And those are equal! It doesn’t matter which one you choose, but I find Judd and Kenny to be a little more straightforward.


library(mediation)

medTest <- mediate(pathA, pathFull, sims = 50, 
                    boot = TRUE, 
                    treat = "NARCOTICS_ARR_CNT", 
                    mediator = "WEAPONS_ARR_CNT")
summary(medTest)

Causal Mediation Analysis 

Nonparametric Bootstrap Confidence Intervals with the Percentile Method

               Estimate 95% CI Lower 95% CI Upper p-value    
ACME            -0.1206      -0.2246        -0.05  <2e-16 ***
ADE              2.8072       2.2102         3.18  <2e-16 ***
Total Effect     2.6866       2.0442         3.09  <2e-16 ***
Prop. Mediated  -0.0449      -0.0926        -0.02  <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Sample Size Used: 9040 


Simulations: 50 

Some of those values should look pretty similar to what we have already done! ACME (Average Causal Mediation Effect) is the indirect effect where we multiplied or substracted our coefficients. ADE (Average Direct Effect) is from our full model. The total effect is what we found from pathC. We are most focused on ACME here, and this significant effect would indicate that we have a significant indirect effect.

We can learn a lot through mediation (it can also be combined with interactions to form mediated moderation or moderated mediation), but it should be used with a great deal of though; it is not just something you start testing models to see how they work. Typically, you need to have a pretty good theoretical reason to dive into mediation.

Effect Sizes

Effect sizes, in conjunction with our p-values, will provide a really good idea about the strength of the difference.

With regard to effect sizes, you will most commonly come across Cohen’s d – it is generally used for t-tests.

Computationally, it is pretty simple:

\[ \frac{\mu_1 - \mu_2}{\sigma}\]

There is also an expanded version:

\[ d = \frac{\mu_1-\mu_2}{\sigma_{pooled}} \]

\[ SD_{pooled} = \sqrt{\frac{\sigma_1^2 + \sigma_2^2}{2}} \]

We are subtracting the mean of one group from another and dividing by the standard deviation.


library(dplyr)

crimeScoreGender[, .(mean = mean(SSL_SCORE), 
                     sd = sd(SSL_SCORE), 
                     n = .N), 
                 by = SEX_CODE_CD]

   SEX_CODE_CD     mean       sd      n
1:           M 278.6891 59.52397 302320
2:           F 283.4600 52.74889  96307

# crimeScoreGender %>%
#   group_by(SEX_CODE_CD) %>%
#   summarize(mean = mean(SSL_SCORE),
#             sd = sd(SSL_SCORE),
#             n = n())

sd(crimeScoreGender$SSL_SCORE)

[1] 57.99564

We can do it by hand:


(283.46-278.689) / 57.99564

[1] 0.0822648

And with the pooled method:


sdPooled <- sqrt((52.74889^2 + 59.52397^2) / 2)

(283.46-278.689) / sdPooled

[1] 0.08483505

Or use things already built:


library(compute.es)

tes(t = 23.674, n.1 = 96307, n.2 = 302320)

Mean Differences ES: 
 
 d [ 95 %CI] = 0.09 [ 0.08 , 0.09 ] 
  var(d) = 0 
  p-value(d) = 0 
  U3(d) = 53.49 % 
  CLES(d) = 52.47 % 
  Cliff's Delta = 0.05 
 
 g [ 95 %CI] = 0.09 [ 0.08 , 0.09 ] 
  var(g) = 0 
  p-value(g) = 0 
  U3(g) = 53.49 % 
  CLES(g) = 52.47 % 
 
 Correlation ES: 
 
 r [ 95 %CI] = 0.04 [ 0.03 , 0.04 ] 
  var(r) = 0 
  p-value(r) = 0 
 
 z [ 95 %CI] = 0.04 [ 0.03 , 0.04 ] 
  var(z) = 0 
  p-value(z) = 0 
 
 Odds Ratio ES: 
 
 OR [ 95 %CI] = 1.17 [ 1.16 , 1.19 ] 
  p-value(OR) = 0 
 
 Log OR [ 95 %CI] = 0.16 [ 0.15 , 0.17 ] 
  var(lOR) = 0 
  p-value(Log OR) = 0 
 
 Other: 
 
 NNT = 39.34 
 Total N = 398627

mes(m.1 = 283.46, m.2 = 278.689,
    sd.1 = 52.74889, sd.2 = 59.52397,
    n.1 = 96307, n.2 = 302320)

Mean Differences ES: 
 
 d [ 95 %CI] = 0.08 [ 0.08 , 0.09 ] 
  var(d) = 0 
  p-value(d) = 0 
  U3(d) = 53.28 % 
  CLES(d) = 52.32 % 
  Cliff's Delta = 0.05 
 
 g [ 95 %CI] = 0.08 [ 0.08 , 0.09 ] 
  var(g) = 0 
  p-value(g) = 0 
  U3(g) = 53.28 % 
  CLES(g) = 52.32 % 
 
 Correlation ES: 
 
 r [ 95 %CI] = 0.04 [ 0.03 , 0.04 ] 
  var(r) = 0 
  p-value(r) = 0 
 
 z [ 95 %CI] = 0.04 [ 0.03 , 0.04 ] 
  var(z) = 0 
  p-value(z) = 0 
 
 Odds Ratio ES: 
 
 OR [ 95 %CI] = 1.16 [ 1.15 , 1.18 ] 
  p-value(OR) = 0 
 
 Log OR [ 95 %CI] = 0.15 [ 0.14 , 0.16 ] 
  var(lOR) = 0 
  p-value(Log OR) = 0 
 
 Other: 
 
 NNT = 41.96 
 Total N = 398627

Power Analysis

Rules of thumb have been around for a long time and have changed over the years – maybe you learned that you needed 20 rows per predictor, or maybe even 50 rows per predictor. Instead of trusting outdated advice, use actual science to determine how many people you need to find if a difference exists.

We need three of the following parameters:

We should always be doing this a priori.

Power

Power is ability to detect an effect. In NHST words, we are trying to determine if we correctly reject the null hypothesis.

Putting It All Together

Let’s use the pwr package.


library(pwr)

pwr.f2.test(u = 1, v = NULL, f2 = .05, power = .8)

     Multiple regression power calculation 

              u = 1
              v = 156.9209
             f2 = 0.05
      sig.level = 0.05
          power = 0.8

In the function:

Power is typically set at .8, because it represents a 4 to 1 trade between Type II and Type I errors.

Different Test, Different Power Tests

We just did a test for a linear regression model.

Here is one for a t-test:


tPower <- pwr.t.test(n = NULL, d = 0.1, power = 0.8,
                    type = "two.sample", alternative = "greater")

plot(tPower)

Let theory be your guide (but be realistic).